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ABSTRACT 

We describe a simple test of the spatial uniformity of an ensemble of discrete events. 
Given an estimate for the point source luminosity function and an instrumental point 
spread function (PSF), a robust upper bound on the fractional point source contri- 
bution to a diffuse signal can be found. We verify with Monte Carlo tests that the 
statistic has advantages over the two-point correlation function for this purpose, and 
derive analytic estimates of the statistic's mean and variance as a function of the point 
source contribution. As a case study, we apply this statistic to recent gamma-ray data 
from the Fermi Large Area Telescope (LAT), and demonstrate that at energies above 
10 GeV, the contribution of unresolved point sources to the diffuse emission is small 
in the region relevant for study of the WMAP Haze. 

Key words: methods: statistical - gamma-rays: diffuse background. 



1 INTRODUCTION 

Statistical tests of isotropy have a long history in astronomy. 
A common question is "What fraction of the observed emis- 
sion could originate from unresolved point sources?" For ex- 
ample, possible point source contributions t o the extragalac - 
tic X-ray background were investigated bv lScheuerl (| 19741 ) . 
and the small-angle power spectrum of the cosmic far infra- 
red back ground has been used to estimate the isotropic com- 
ponent (|Kashlinskv et all Il996t ). More recently the Auger 
team tested the isotropy of ultra-high energy cosmic ray 
events by cross correl ating with posit i ons of know n active 
galactic nuclei (AGN; [Abraham et ail l|2007l . 120081 )) to pro- 
vide information on their origin. However in the absence 
of an appropriate external catalog, such cross-correlation 
methods cannot be used, motivating consideration of a more 
general approach. 

In some cases a detector provides binned counts (e.g. 
pixels in a CCD); in other cases, photon event directions are 
reconstructed in some other way (e.g. a gamma-ray pair con- 
version telescope). In the latter case, it is desirable to apply 
statistics that do not require binning of the data, as binning 
introduces additional arbitrary parameters into the prob- 
lem. In the limit of low flux density, where the mean density 
of photon events (hereafter, "events") is much less than one 
per PSF, explicit detection of point sources may become im- 
practical, and estimation of the unresolved point source flux 
becomes especially difficult. In some cases, the two-point 
corre lation function, or some modified form (e.g. lAve et all 
2009), is used as a test. However, the Fourier transform of a 
field of point sources has significant phase correlation, and a 
two-point function (or a power spectrum) discards this phase 



information. Higher order correlation statistics capture it, 
but are somewhat cumbersome to use. In the following, we 
describe a statistic that is easy to understand and evaluate, 
and is optimized to address this question, particularly in the 
case of fairly sparse data sets with (on average) < 1 event 
per PSF circle. 

The key insight is that if a substantial fraction of the 
photons come from point sources, it is much more likely that 
two photons appear within one PSF of each other than in the 
diffuse case. This is true even if the expected integrated flux 
is of order one count. Likewise, the number of PSF circles 
containing no counts is larger if point sources contribute. 
These considerations motivate us to define a ratio between 
the fraction of "isolated events" and the fraction of "empty 
circles." This ratio is very closely related to the fraction of 
diffuse emission, and can be calibrated with Monte Carlo 
simulations for specific choices of instrumental parameters 
and a putative luminosity function. The two-point function, 
in contrast, is weighted by the density squared and is not 
proportional to the desired quantity. 

In the following sections we define the statistic, estimate 
its variance, show how it behaves in various limits, generalize 
it to the case where many events appear in every PSF circle, 
and show a practical application to recently released data 
from the Fermi Gamma-ray Space Telescope. 



2 DEFINITION OF THE STATISTIC 

For each event, we consider the number of neighbouring 
events within some test radius r; the natural choice of r is 
determined by the PSF of the detector. Let the fraction of 
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Figure 1. In each panel, photon events (dots) arc cither isolated (solid blue circles) or not (dashed blue circles). Random circles are 
either empty (solid red) or not(dashed red). In each case, 100 events and 100 random circles are shown, so the ratio, R, can be visualized 
here as the number of solid blue circles divided by the number of solid red circles. In practise, one uses a large number of random circles 
to reduce noise. The panels contain either no point sources (upper left), or 15% (upper right), 50% (lower left), or 90% (lower right) 
point source flux in sources located at (0.3,0.3) and (0.7,0.7). 



events with zero neighbours ("isolated events") be denoted 
ni . Now consider a random distribution of points within the 
signal region, and consider the number of events within r of 
each of these points. Let the fraction of points with zero 
events within radius r ("empty circles") be denoted tie (the 
number of randomly distributed points should be as large 
as is computationally feasible, to reduce Poisson error in 
the fraction tig). 

Fig. [T] illustrates this idea. If the signal photons are 
randomly distributed, as in the case of a uniform diffuse 



signal, ni ~ He (up to Poisson fluctuations). If the events 
are clustered on the scale r, on the other hand, ni falls (as 
more events have close neighbours) and tie rises (as the 
events are clumped, more of the signal region contains no 
events at all). 



If a large fraction of the photon counts originate from 
unresolved point sources, we expect significant correlations 
between photon positions, even if the expected counts from 
a point source are < 1. Consequently, the ratio R — ni /tie 
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serves as a simple measure of the fraction of flux originating 
from point sources as opposed to uniform diffuse emission. 

Larger-scale non-uniformities, such as a gradient in the 
distribution of diffuse photons, have only a small effect on 
R, as we show explicitly with Monte Carlo tests in fj6] In re- 
gions of high event density both ni and tie are suppressed, 
so the suppression largely cancels out in the ratio. However, 
such a gradient will tend to slightly lower R in the diffuse 
limit, since the points relevant to computing tie are uni- 
formly distributed, whereas those relevant to computing nj 
are concentrated in regions of high flux (this in turn can re- 
duce the gradient of R with respect to the fraction of diffuse 
flux, since in the point-source-dominated limit 7? is small 
but independent of any inhomogeneity in the diffuse emis- 
sion). The point source emission derived from R is thus best 
interpreted as an upper bound, in cases where the diffuse 
emission is suspected to have significant spatial variation. 



3 ANALYTIC ESTIMATE FOR THE MEAN 
AND VARIANCE 

To understand the behaviour of the ratio J? as a function of 
the diffuse emission fraction, consider a related but simpler 
problem, where we treat the signal region as a grid and count 
the number of events in each cell. In this case an "isolated" 
event is one with no other events in the same cell, and R is 
simply the ratio of the fraction of events which are isolated 
to the fraction of cells which are empty. The total number 
of cells is N, n of which are empty and m of which have a 
single event. 

The probability of any given cell being empty is po, 
and of containing a single event, pi. That is, (n) — poN, 
and (m) = piN. The joint probability of having precisely n 
empty cells and m single-event cells is 

N\ 



(m,n) = p Pi (1-Po-pi) 



n\m\{N — m — n)\ 



(1) 



The fraction of isolated events is (ni) = m/N cvcnt , while 
the empty fraction is {tie) = n/N. In this case R = m/(n + 
l)x TV/iVevcnt. As we shall see, (m/(n + I)) is an unbiased 
estimator of the probability ratio pi/po, which corresponds 
to R as defined above in the limit of large n, and is equivalent 
to the fraction of the flux originating from diffuse emission 
if the point sources are sufficiently bright. 

The expectation value of m/(n + 1) is given by, 
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Consequently, neglecting terms suppressed by a large power 
of (1 — po), the mean fi and variance a 2 of m/(n + 1) in the 
case of large N are given by, 



A 4 — Pi/ 



Pl (l+ Pl /p )/Npl 



(4) 



The corresponding quantities for R are obtained by a rescal- 
ing by N/N cvcnt and (N/N event ) 2 , respectively: 



(R) = 



Pi N ( <tr\ 2 _ 1 +P1/P0 _j_ 

PO jVcvcnt ' \{R)J ~ JVpi ~ m" U 



It is unsurprising that the fractional uncertainty in R is ap- 
proximately 1 / y/m where m is the number of isolated events. 

In the limit where all the flux is diffuse with mean rate 
A, then pi = Ae~ A , po = e _A , and {N cvcn t} = iVA, so R 
has mean 1 and variance (1 + A)e A /7VA. Now consider the 
addition of point sources. We can approximate the effect of 
adding point sources by choosing T cells which each gain > 1 
events. Let the total number of added events be T. Then on 
average, a fraction T/N of the isolated events will no longer 
be isolated, and a fraction T/N of the empty cells will no 
longer be empty (by assuming each cell that gains events 
gains > 1, we ensure that empty cells will not gain isolated 
events). 

Thus when point sources are added, both po and pi are 
multiplied by (1 - T/N), and Advent becomes 7VA + V. Thus 
the mean of R becomes NX/ (NX + F) = diffuse flux / total 
flux, and the variance becomes, 



a 2 (R) 



(l + A)e A 
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(6) 



diffuse flux/total flux 



fraction of pixels with no point sources 



So in this simple model, the mean of the ratio R is pre- 
cisely the diffuse flux divided by the total flux. The variance 
in R is generally Poisson, scaling as l/(total photon counts) 
provided the number of counts per pixel from diffuse emis- 
sion is < 1, but is reduced when the diffuse flux is much 
smaller than the total flux. (It is possible for the diffuse flux 
to be much less than the total flux even when most cells 
do not contain point sources, but not the converse, unless 
A is rather large, in which case the sample size or the cell 
size should be reduced.) The variance grows rapidly as A 
becomes greater than 1: if the cell size (corresponding to 
angular resolution) is too large relative to the amount of 
diffuse emission, R is not a good measure of the fraction of 
diffuse flux (however, see £|3.1[) . 

If there is a significant contribution from weak point 
sources that may only add a single isolated event to an 
empty cell (or alternatively, if the cells are sufficiently small 
that this is common even for stronger point sources), then 
the mean of R is no longer given simply by the fraction of 
diffuse flux. For example, suppose S cells gain exactly one 
event when point sources are added, and T cells in total are 
affected by the addition of point sources (so T > S): then 
the addition of point sources sends po — > (1 — T/N)po and 
pi ->• (1 - T/N) Pl + (S/N)p . Then we obtain, 
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«rm = ( 1 ^ eA ( 1 + A + ^) ( 1 + ^t)) 
a[ ' ^iVA + r^ (1 + r/NX) (i — t/n) 

(8) 

If the number of cells containing single 1-photon point 
sources is small compared to the number of cells unaffected 
by point sources (S <C N — T), and also the number of pho- 
tons from these weak point sources is small compared to 
the number of diffuse events from cells unaffected by point 
sources (S <C X(N — T)), we recover the previous result. If 
the second condition fails to hold (as occurs in the limit of 
low diffuse emission, independent of the point source lumi- 
nosity function), then (R) asymptotes to (flux from isolated 
1-photon point sources) /(total flux) xl/(l— T/N), as A — > 
with the total flux held constant. If this limiting value is ^ 1, 
then R has no discriminatory power, and the situation is not 
improved by higher statistics: this is simply the statement 
that there is no difference between diffuse emission and a 
very large number of very faint uniformly distributed point 
sources. 



3.1 Extension to the case of large A 

One region of parameter space in which this test breaks 
down is where A > 1. However, a simple generalization of 
the statistic can be useful in this case. Suppose we make a 
histogram of the number of nearest neighbours each event 
possesses, and measure the peak of the histogram to be some 
number of neighbours n cr i t . Let us redefine nj as the frac- 
tion of events with n cr i t or fewer neighbours, and tie as the 
fraction of points with n cr j t or fewer events within the test 
radius r. Then for Monte Carlo realizations of diffuse flux 
plus some randomly distributed point sources, we can again 
measure the ratio Ui/ue and use it as a measure of how 
correlated the photon events are (see i]6.2[) . 



4 THE POINT SOURCE LUMINOSITY 
FUNCTION 

Clearly, the sensitivity of the test depends critically on the 
fraction of point sources which contribute at most one count 
to the data, which is determined by the point source lumi- 
nosity function. In the limit where the point source flux is 
dominated by (a very large number of) uniformly distributed 
sources which each produce an average number of counts 
<C 1, unresolved point sources are practically indistinguish- 
able from diffuse emission. 

In the Monte Carlo tests which follow, we treat the lu- 
minosity function as some unbroken power law dN/dS = 
S~ a between integrated flux limits S m i n and S m ax- S m ax 
is bounded above by the faintest point sources which can 
be resolved and masked out. S m m, on the other hand, is 
not determined by known properties of the experiment, and 
limits derived with a particular value of S m in should be in- 
terpreted as placing a limit on the contribution from point 
sources with average luminosity above S m in. 

The effect of S m in on the behaviour of R depends on 
the value of the spectral index a: if a > 2, then most of the 
flux originates from the faintest point sources, and the total 
flux diverges at low luminosity. In this case, the power law 



generally breaks to a much shallower slope at some low lumi- 
nosity, meaning that most of the flux originates from point 
sources with luminosities close to the break: we can approx- 
imate this behaviour simply by cutting off the luminosity 
function at S mm . If a < 2 then most of the flux is concen- 
trated in the brighter point sources and changing S'min has 
little effect on R. 

To estimate the relevant range of a we examine stud- 
ies of known populations of point sources. The gamma- 
ray luminosity functions of AGN contributing to the ex- 
tragalactic gamma ray background have been studied us- 
ing the first three months of data from the Fermi LAT 
l)Abdo et al.ll2009aT ) . The luminosity function for BL Lac ob- 
jects was found to be well described by a single power law 
with a — 2.17 ± 0.05. The luminosity function for flat spec- 
trum radio quasars (FSRQs) was well described by a power 
law with a = 2.58 ± 0.19 at high redshifts (z > 1), indicat- 
ing that at high redshifts the Fermi LAT is sampling the 
bright (steep) end of the luminosity function, but at z < 1 
the best-fit value of the slope was a — 1.56±0.10. The X-ray 
luminosity functions of the same classes of objects have been 
studied bv lPadovani et all l)2007h . with BL Lac objects mea- 
sured to have a = 2.12 ± 0.16, and a = 1.6 - 1.9 for FSRQs 
at z < 1. 

The X-ray luminosity fun ction of high-mass X -ray bi- 
naries (HMXBs) was found bv lGrimm etail l|2003l ) to have 
slope a — 1.61 ± 0.12. St udies of low-mass X-ray binaries 
(LMX Bs) in Centaurus A (jVoss fc Gilfanovll2006l : IVoss et all 
2009) yield a slope of a = 1.8 — 2.0 at high luminosity, 
flattening to a ~ 1.2 at lo w luminosity, in agreement with 
earlie r studies of LMXBs ijKim fc Fabbianolbooi iGilfanovl 
l2004h . 

In the examples in the following section, therefore, we 
take our two benchmark models to have a = 1.8, 2.2, and 
also demonstrate the effect of varying a between 1.5 and 
3.0. Smaller values for a improve the ability of our statistic 
to distinguish point sources from diffuse emission, simply 
because more of the flux originates from bright sources. 



5 MONTE CARLO EXAMPLES 

To examine the usefulness of the ratio R as a measure of the 
point source flux, we employ a Monte Carlo approach. We 
consider a given angular "signal region" , and a smaller "sig- 
nal window" within that region (the purpose of this distinc- 
tion is to eliminate edge effects) . Within the signal region we 
generate a uniform random distribution of point source lo- 
cations, for a certain number of sources (which we vary, as a 
proxy for varying the total flux from point sources) . The ex- 
pected number of counts for each point source is drawn from 
a power law distribution with spectral index a, with cutoffs 
at a minimum expected number of counts S'min, and a max- 
imum expected number of counts S max . For our benchmark 
models we take (1) a "pessimistic" set of parameters a = 2.2, 
S m ax = 10, Smin = 0.1, and (2) an "optimistic" set of pa- 
rameters a = 1.8, S m ax = 100, Smin = 1 (see Tabled]). Once 
the expected counts from each source have been obtained, 
the number of counts actually observed from each source 
is determined by a Poisson draw; their angular distribution 
is determined by the detector PSF, which for simplicity we 
model as a Gaussian with standard deviation a. 
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We compute the total number of counts lying within 
the signal window and originating from the point sources, 
and subtract this quantity from the desired total number of 
counts in the signal window. The result approximates the 
number of diffuse events we wish to generate in the signal 
window. We then generate a uniform random distribution 
of diffuse events in the signal region, with the total number 
of events given by the desired diffuse counts in the signal 
window rescaled to the greater area of the signal region. 
The total photon distribution in the signal window (diffuse 
+ point sources) then has approximately the correct number 
of events. We can also incorporate non-uniformities in the 
diffuse flux at this point. 

Having produced our test data, we need to compute the 
ratio R. For each event within the signal window, we use 
the publicly available IDL routine spherematchQ to find its 
neighbours in the set of events in the larger signal region. 
This eliminates edge effects, i.e. a spuriously high number of 
isolated events at the edges of the signal window. We repeat 
the process for a random distribution of points within the 
signal window, to compute tie, and plot the resulting R 
against the fraction of flux in the signal window originating 
from point sources. 

For these examples, we take the default signal window 
to be the region \l\ < 15, \b\ < 15, the PSF to be 0.2°, and 
the total number of counts to be 3000, corresponding to a 
mean A = 0.4 events per PSF circle. Note that the statistic 
is insensitive to the shape of the signal window, and if we 
rescale the PSF by some factor a and the area of the signal 
window by a 2 , and hold the number of events constant, then 
this is just equivalent to a unit redefinition and does not 
change the results. 

When we refer to the "PSF", we mean the standard 
error la. Optical astronomers often use the full width at 
half max (FWHM), which for a Gaussian equals 2.355cr and 
contains 50% of the flux. Gamma-ray astronomers often use 
the radius of 68% or 95% containment (1.51a or 2.45<r, re- 
spectively). The R statistic is not greatly sensitive to the 
exact choice of r. We will demonstrate the effect of varying 
the PSF while holding the other parameters constant (thus 
changing A, the mean number of events per PSF circle), and 
of increasing the number of counts while holding A constant, 
by increasing the signal window. We will also generally as- 
sume that the test radius r is equal to la, but show the 
effect of using a different test radius, for both benchmark 
luminosity functions. 

In the plots that follow, we will describe the "sensitiv- 
ity" of this test by two representative measures: 

(i) The maximum value of the point source flux fraction 
consistent with R = 1, within the 95% confidence limits. 

(ii) Bounds on the point source flux fraction when the 
true point source contribution is half the signal, obtained 
by averaging the limits which would be obtained from an 
R-measurement over the histogram of J?-values produced in 
this scenario. 

The statistic is most powerful where the bound in (i) is small 



The IDL routines used in this analysis are avail- 
able as part of v5_4_8 of the IDLUTILS product at 
http : //sdss3data. lbl . gov/sof tware/idlutils 



Table 1. Benchmark parameters for luminosity functions 
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and the limits in (ii) are close together. If either the variance 
of R (for fixed point source flux) becomes large, or the mean 
of R is slowly varying with respect to the point source flux 
fraction, then both of these measures will blow up. We also 
overplot the sensitivity estimates obtained from the analytic 
approximation for R in the grid model described in Sj3] using 
the relations in Appendix [A] for comparison to the results 
of the MC runs. 



6 RESULTS 

We find a strong linear relationship between R and the frac- 
tion of the flux due to diffuse emission, as shown in Fig. 
for the two benchmark sets of parameters. As expected, the 
"optimistic" benchmark parameters render R more sensitive 
to the fraction of diffuse emission. 

6.1 Dependence on the point source luminosity 
function and the effect of improved statistics 

Fig. [3] shows the effect of varying the spectral index of the 
point source luminosity function, for the two benchmark 
choices of S'min, S'max- As expected, smaller values of a yield 
better performance for R as a precise estimator of the diffuse 
flux. For the pessimistic benchmark model, when a ~ 3.0 the 
method has lost most of its discriminatory power: however, 
for more well-motivated power laws a < 2.5, useful limits 
can still be obtained. For the optimistic benchmark model, 
R can be used to place strong limits on the diffuse emis- 
sion even for very steep power laws a ~ 3.0, largely because 
the assumed break in the power law at S m m = 1 avoids the 
scenario of emission dominated by many very faint point 
sources (S <C 1). 

As a check on our understanding of the method, we can 
increase the area of the signal window and also the num- 
ber of counts, holding the density of events constant. This 
change would be expected to have no effect other than re- 
ducing the variance in R, since a(R) is expected to scale as 
1/%/number of counts. Fig.[3]demonstrates the improvement 
in sensitivity. 

6.2 Dependence on the PSF and test radius r 

The choice of the test radius r significantly impacts the qual- 
ity of the results. If r is increased sufficiently that the ex- 
pected number of neighbours for each event is ^> 1, then 
the number of isolated events and empty circles both be- 
come small in the diffuse limit, and the Poisson fluctuations 
in R become very large. In the limit where r <C PSF, even 
events in bright point sources may qualify as "isolated" and 
R loses its power to discriminate between point sources and 
diffuse signal. In the large-r case, a better result with much 
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Figure 2. The isotropy ratio R, in Monte Carlo simulated data, in (left) benchmark model 1 and (right) benchmark model 2. See Table 
[T]for definitions of the benchmarks. Lines are 5% and 95% confidence bounds. 
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Figure 3. The sensitivity of R as a function of the spectral index of the point source luminosity function, a (far left: benchmark model 
1, center left: benchmark model 2), and as a function of the signal area while the density of counts is held constant (center right: 
benchmark model 1, far right: benchmark model 2). As described in the text, solid lines (bounding the shaded area) indicate the average 
90 % confidence bounds on the point source flux fraction from this test when the true fraction is 0.5; the dashed line indicates the 
95% confidence upper limit on the point source flux fraction where R = 1. Dotted lines indicate the analogous results for the analytic 
approximate calculation. 



less noise can be obtained by the generalized method dis- 
cussed in H3.ll These effects are displayed in Fig. 2] 

The range of r in which R provides a precise estimate 
of the point source flux is given approximately by, 



PSF, the breakdown in the method at large PSF, and the 
performance of the generalized method. 



PSF< r< 



area of window 
7T x number of events 



(9) 



If this range is large, r can be varied substantially without 
much adverse impact on the performance of R as a measure 
of diffuse flux. For the parameters employed here, the per- 
mitted range of r is quite narrow and the effects discussed 
above are pronounced. 

Where this allowed range vanishes, as the average num- 
ber of counts per PSF circle exceeds 1, this method breaks 
down as discussed in Sj3] and we should instead employ the 
generalization discussed in cj3.ll Fig. [4] shows the slight im- 
provement in the sensitivity of R in the case of a smaller 



6.3 Inhomogeneity in the diffuse flux 

The "diffuse" part of the sample may vary spatially, either 
because of true spatial variation of the signal, or a non- 
uniform instrumental sensitivity or exposure. We show that 
R is insensitive to such variations by introducing a tilt in 
the diffuse photon distribution of up to a factor of 20 (that 
is, the density of events is 20 times lower at one edge of the 
signal region than at the other) and observing that it has 
little effect on R, as shown in Fig. [5] For large tilts, there is 
a noticeable downward shift in R in the diffuse limit. 
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Figure 4. The sensitivity of Ft as a function of the test radius r, holding the PSF lcr constant at 0.2° (far left: benchmark model 1, 
center left: benchmark model 2), and as a function of the PSF size, fixing the test radius r to be equal to lcr for the PSF (center right: 
benchmark model 1, far right: benchmark model 2). The top row uses the standard form of the statistic whereas the bottom row uses 
the modified form ([[3Tj to prevent the catastrophic failure at large r. As described in the text, solid lines (bounding the shaded area) 
indicate the average 90 % confidence bounds on the point source flux fraction from this test when the true fraction is 0.5; the dashed 
line indicates the 95% confidence upper limit on the point source flux fraction where Ft = 1. In the top row, dotted lines indicate the 
analogous results for the analytic approximate calculation; note that the truncation approximation made to derive the analytic result 
breaks down for r -C lcr PSF. 
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Figure 5. The effect of adding a 6-dependent tilt to the diffuse photon flux, in Monte Carlo simulated data, for (far left) benchmark 
model 1, tilt factor 5, (center left) benchmark model 1, tilt factor 20, (center right) benchmark model 2, tilt factor 5 and (far right) 
benchmark model 2, tilt factor 20. 



7 COMPARISON TO THE TWO-POINT 
FUNCTION 

In essence, R is a simple measure of the angular correlations 
between event positions, so it is reasonable to ask how it 
differs from the two-point correlation function. Consider the 
case of a single bright point source, compared to two (well 
separated) point sources with half the luminosity. The two- 
point correlation functions for these two situations are quite 
different, although the total flux from the point sources is 
the same, because the number of pairs in a given source 
scales as flux squared. The estimator R, on the other hand, 
to a first approximation does not probe correlations inside 
the point sources, and is dependent only on the fraction of 
diffuse flux, rather than the details of the sources. 



This effect tends to reduce the variance of R, relative 
to the two-point function, as the flux from point sources in- 
creases and the luminosity function becomes shallower (i.e. 
brighter point sources, with (no. of counts) 2 ^> no. of counts, 
contribute more of the signal); the R statistic is also less sen- 
sitive to the luminosity function parameters than the two- 
point function. In the limit where the annulus width is large 
and almost all the non-isolated events have only a single 
neighbor, i.e. they form a single pair in the calculation of 
the two-point function (which is the case when the luminos- 
ity function is steep and the point source flux is dominated 
by faint sources, or when there are simply very few point 
sources), the two statistics capture essentially the same in- 
formation and their performance is very similar. 
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We can directly compare these two statistics as mea- 
sures of the point source flux. We employ the unbi- 
ased estimator f or th e two-point function described by 
lLandv fe Szalavl (|l993h . and compute its value as a func- 
tion of the point source flux contribution (we also scan over 
the annulus width A9), within the Monte Carlo framework 
described previously. Applying the sensitivity measures pre- 
viously described, the results for the two benchmark param- 
eter sets are shown in Fig. [6] As expected from the discus- 
sion above, we see that in the benchmark 1 case, the results 
are very similar to those previously found, whereas in the 
more "optimistic" benchmark 2 case, the results are simi- 
lar in the diffuse limit, but the R test provides much better 
bounds when the true point source flux fraction is 0.5. 



In all cases we take the test radius r to be equal to the PSF, 
and generate 10 6 random points to determine tie- 

Fig. [7] shows the 95% confidence limits on the fraction 
of flux originating from point sources, as a function of Smin 
and the luminosity function spectral index a. We see that 
even with fairly pessimistic assumptions for the luminosity 
function, as in the "Benchmark I" case, the measured values 
of R in the Fermi data are outside the 95% confidence limits 
if the flux from sources with > 0.1 (expected) counts /year 
exceeds ~ 15% of the total; the limits can be significantly 
stronger if the assumed point source luminosity function is 
shallower or extends to higher flux. Note also that we have 
not subtracted resolved point sources in this region; there 
are no known point sources in this region of the sky in the 
Fermi three- month bright source list q 



8 AN EXAMPLE APPLICATION: HIGH 
ENERGY Fermi LAT DATA FROM THE 
HAZE REGION 

The Fermi Gamma-ray Space Telescope has recently re- 
leased all-sky photon data from its first year of operation 
0. The diffuse emission measured by Fermi may include sig- 
natures of new physics, such as dark matter annihilation or 
decay, or photons from new classes of astrophysical sources. 
In investigating the origin of the diffuse gamma rays, it will 
be necessary to estimate what fraction of the observed flux 
could be due to unresolved point sources. The statistic pre- 
sented here is well suited to address this question, especially 
at high energy where the PSF of the Fermi LAT is small 
and the count rate is low. 

A microwave excess termed the "WMAP Haze" has 
been observed in the inner 25° of the Galaxy, and at- 
tributed to synchrotron radiation from some new pop- 
ulation of 10-1000 GeV electrons (and/or positrons) 
jPobler fc Finkbeinerll2008l '). Recent cosmic ray experiments 
have also measured a rise in the positron fraction at 10—100 
GeV (|Adriani et al.ll2009h . and a harde ning in the e + e 
spectrum at 300 GeV - several TeV l|Chang et all 120081 ; 
lAbdo et all 12009^. consistent with a new source of high 
energy e + e~ . The Fermi LAT can constrain any such new 
source of electrons by searching for gamma rays from inverse 
Compton scattering of the electrons on starlight. 

As an example of how this statistic can be applied, we 
consider the Class 3 (diffuse class) events measured by the 
Fermi LAT with energies between 10-100 GeV, in the re- 
gion of the sky optimized for study of the WMAP Haze, 
defined in Galactic coordinates by |Z| < 15, —30 < b < —10. 
There are 1146 such events in this signal window. We take 
the signal region \l\ < 18, -40 < b < -8. Above 10 
GeV the 68% containment radius of the LA T is < 0.2° 
ijRando fc the Fermi LAT Collaboration! [20091 ). which cor- 
responds to a PSF of ~ 0.13° in the sense that we have 
used (ID a for the Gaussian distribution of photons from a 
single point source). We consider three values of the PSF: 
the best-estimate upper bound of 0.13°, and also 0.1° and 
0.2°, to demonstrate the possible effect of uncertainties in 
the PSF, or a varying PSF over the energy range of interest. 



9 CONCLUSION 

We have introduced a simple and easily calculable statis- 
tic that linearly traces the fraction of flux arising from dif- 
fuse emission, as opposed to unresolved point sources. The 
statistic is quite insensitive to even pronounced large-scale 
anisotropies in the diffuse emission, such as might originate 
from the proximity of a bright region or angular variation in 
the detector exposure. The linear response of this statistic 
to flux originating from point sources, and its smaller vari- 
ance, make it superior to the two-point correlation function 
as a tracer of emission from unresolved point sources. 

The sensitivity of the statistic to point source emission 
naturally depends on the luminosity function of the point 
sources, as a sufficiently steep power law extending to suf- 
ficiently small luminosities is strictly indistinguishable from 
diffuse emission. However, the statistic retains discrimina- 
tory power for spectral indices up to a ~ 3, with a low- 
luminosity cutoff corresponding to an average of 0.1 counts, 
and assuming all point sources with average luminosity > 10 
counts are resolved and removed. Known luminosity func- 
tions for astrophysical point sources generically have shal- 
lower slopes than this limit. 

When the average number of events per PSF circle ex- 
ceeds 1, the original form of the statistic breaks down: how- 
ever, we have described a simple generalization suitable for 
this case, and demonstrated its efficacy. Increasing the num- 
ber of counts by taking additional sky regions into account 
(i.e. without increasing the density of events) improves the 
variance by the usual 1/N even t Poisson factor. This statistic 
generalizes readily to higher dimensions; possible applica- 
tions include the study of void statistics (|FrvHl986l ). 

As an example, we have applied this statistic to Class 3 
(diffuse class) photon data from the Fermi LAT in the angu- 
lar region relevant for study of the WMAP Haze, at energies 
of 10-100 GeV. We find that even with rather pessimistic as- 
sumptions for the point source luminosity function, at most 
~ 15% of the emission in this region can be attributed to 
unresolved point sources with average luminosities of 0.1+ 
counts / year, and the results are consistent with 100% dif- 
fuse emission. 

We wish to acknowledge helpful conversations with 



See http : //f ermi .gsf c .nasa. gov/ssc/data/ 3 http: / /f ermi .gsf c.nasa. gov/ ssc/data/access/lat/bright_src_list/ 
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Figure 6. Monte Carlo realizations of the two-point correlation function for the angular region and number of events described in [[5] 
varying the annulus width as a fraction of the test radius (the fractional annulus width is denoted A8) and taking r=0.2°=lcr of the 
PSF, for (left) benchmark model 1, and (center) benchmark model 2. The right-hand panel shows an example of the mean and 5% and 
95 % quantiles for the two-point function as a function of the point source flux fraction, in benchmark model 2 with Ad = 0.5. 
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Figure 7. Bounds on the point source flux fraction in Fermi LAT data in the "Haze" region at 10-100 GeV, from the isotropy ratio R 
in Monte Carlo simulated data, as a function of the spectral index of the luminosity function. In all cases 5 ma x = 10. Left: r = 1 a PSF 
= 0.1°, center: r = 1 a PSF = 0.13° (the estimated PSF for back-converting events at 10 GeV), right: r = 1 a PSF = 0.2°. Black solid: 
•Smin = 0.01, blue dotted: 5 m i n = 0.1, red dashed: 5 m i n = 1. 
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© 0000 RAS, MNRAS 000, 000-000 



10 Tracy R. Slatyer and Douglas P. Finkbeiner 



APPENDIX A: THE ANALYTIC GRID MODEL AND THE LUMINOSITY FUNCTION 

The analytic estimates for R and cr(R) derived in ij3]are functions of the fraction of cells which gain exactly 0-1 counts from 
point sources. If we are to compare the analytic estimates to the results of the MC runs, the S and T parameters must be 
expressed in terms of a, S m in and S'max- This is a nontrivial exercise, and in any case the analytic estimates and the MC runs 
should not be expected to agree in detail, since they use different criteria for determining neighboring events (events within 
the same cell vs events within the test radius). The main purpose of the analytic estimates is to demonstrate the general 
scaling behavior of 7? and cr(R), and the regions of parameter space where this test loses discriminatory power. 

Nonetheless, for completeness, we now derive approximate relations between S and T and the parameters of the luminosity 
function. These relations are used to provide an analytic estimate for the sensitivity of R as a measure of the point source 
flux fraction, in Figs. [3J4] Here we only derive the mean values for S and T for a given luminosity function; for any choice of 
the luminosity function, there will be an additional contribution to cr 2 {R) from the variances of S and T, which is not taken 
into account in this analysis. 

Let us first compute the probabilities for a given cell to obtain exactly one or zero counts from the addition of a 

single (randomly placed) point source, denoted and p' 3 respectively. If the point source has an expected contribution of k 
counts, the PSF is assumed to be Gaussian with variance a 2 , the cells are square with dimensions 5 x S, and the point source 
is centered at (xo, yo), then the expected number of counts in a cell with left-hand lower corner (a;, yj) is given by, 

Ay (fc, *<>,»>) = [ 1+S f 3+S e-^ 2 ^ 2 e'^ 2 ^ 2 d X dy, 



x(Evf( Vo - Vj )-Ev{( V °-^-" ) ) ( .VI, 



V2cr J ~ V V2a 

The Poisson probabilities to obtain exactly one and zero counts in this cell, from a point source providing k events, are 
then given by, 

p\>(k,x ,y ) = \ lJ (k,x ,y )e- X ^ xo ' yo \ p^(k,x ,y ) = e - X ^ k ' x °' yo \ (A2) 

Since the Xij's depend only linearly on k, it is straightforward to integrate over the luminosity function: writing Xij = k8, we 
obtain, 



Po(xo,yo) 



is°r k ~ ae ~ k<>dk 

jSmax k - a(Jk 



IS 

na — 1 



(r (1 - a, SminO) - r (1 - a, S max 0)) (1 - a) . 



1 — a ryl — a 



q 1 — a q 

•J max |J, 



mm 



pl J (xo,yo) 



fs27 k ~ akde ~ k6 dk 

jSmax k - adk 

(T- 1 (r (2 - a, S min 8) - T (2 - a,S mMt fl))(l-a) (A4) 



s. 



mill 



Strictly speaking we should now integrate this result with respect to xo and yo; however, this is not analytically tractable. 
Provided a is not too much smaller than 5, it is a good approximation to instead integrate Eq. IA1I over xq and yo within a 
given cell, obtaining an average expected number of counts Xij (k) — k8 for each cell, and then use this result for 6 in Eq. IA4I 
Let us choose our coordinate system so that (xo,yo) lies in a cell with left-hand corner (0,0), and the ij cell has left-hand 
corner (i8,jS). For \i\ or |j| ^> 0, Xij will be negligible, so we can make the further approximation of truncating the sum over 
i,j at some point. In this work we make the approximation that for cells with \i\ or \j\ > 2, p 3 = 1 and p^ = 0. 

We can then write 9 for each cell in terms of three functions of a/8, denoted t, u and v, obtained by averaging 

(1/2) (Erf ( 2 $=f i ) ~ Erf f 3 ^^ )) over x = [0, 5], for i = 0, 1,2 respectively: 
ErfC * (l-e- 62 ^ 2 



(A5) 



In terms of these functions the values of 6 for the relevant cells are: 

i = j = 0, 1 cell, 9 = t 2 
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i = 0, j = ±1 or i = ±1, j = 0, 4 cells, 9 = tu 

i = ±1, j = ±1, 4 cells, 6» = 

i = 0, j = ±2 or i = ±2, j = 0, 4 cells, = tv 

i = ±1, j = ±2 or i = ±2, j = ±1, 8 cells, = uv 

i = ±2, j = ±2, 4 cells, = v 2 . 

Summing over the probabilities p$ (p^ ) for all cells then yields the expected number of cells which gain zero (one) counts 
from the addition of a point source, denoted Eq (Ei). If n sources are added, the probability of any one cell gaining zero 
counts is (Eo/N) n , and the probability of a single count is n(Eo / 7V) n_1 (Ei j TV) (i.e. one source contributes a single count, all 
others contribute zero). Multiplying by the number of cells TV yields the N — T and S parameters respectively. The number 
of sources n is related to the average total emission from point sources, 

51 — OL ol — CK N 

max ^mi n 



• 2 — a 
max 



x mean total counts. (A6) 



In order to compare the analytic result (based on the grid) to the MC results (based on neighbors within a test radius), 
we must also impose a relation between the side length of the grid cells and the test radius. In Figs. 1 3141 we require that the 
area of a grid cell is the same as the area of a PSF circle, i.e. S = v^ir; a different prescription might give better agreement 
between the MC results and the estimates from the grid model. 
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